Plotting Palaeogeographical Maps in R: a Lab Group

Palaeobiology Research Group, School of Earth Sciences, University of Bristol, Bristol, U.K.
benjamin.moon@bristol.ac.uk

1 Introduction

This guide takes you through two different work flows to create palaeogeographical maps and plot occurrences fossil taxa in R. All the code and data are in this GitHub repository: https://github.com/benjaminmoon/palaeomap_example

You can download and run this code directly on your machine, or follow along online at Binder where I’ve created an environment with all of the data you need. Click the link below to open this in your web browser.

The Binder link will open in your browser in the RStudio interface. Open the file palaeomap_labgroup.Rmd from the file browser to get started.

The RStudio interface similar to that you should see in Binder. Open the file palaeomap_labgroup.nb.html to get started.

2 Palaeogeographical Maps

Maps are essential when discussing occurrences of taxa, both in the present and the past. Chuck in the issues associated with continental drift, and it can be difficult to place the wider context of some fossil occurrences. ‘These marine fossils are found at the top of a mountain, but where were the oceans 150 million years ago?’ Here are a couple of recent examples discussing the spread of occurrences, in this case Jurassic marine reptiles.

As a first example, Gao et al. (2019) discuss the occurrences of Jurassic plesiosaurs in China.

Here, Zverkov, Grigoriev, and Danilov (2020) show occurrences of ichthyosaurs in the Early Jurassic.

Reconstructing global palaeogeography has been an ongoing task for many decades, and both local and global reconstructions and compendia have been made (e.g. Bradshaw et al. 1992). But the case of making maps to show many occurrences by hand can be difficult and time-consuming.

A few examples of groups that have made efforts to create palaeogeographical maps at various points through the Phanerozoic include Ron Blakey at Deep Time Maps (https://deeptimemaps.com) and Christopher Scotese leading the Paleomap Project (http://scotese.com).

DeepTimeMaps_Animation_Mollweide_Sample from 42 Degrees North Media on Vimeo.

Palaeogeography of the Late Jurassic, as reconstructed at Palaeomap Project.

While built from models of plate tectonics and continental drift, these are nonetheless individual time points in a continuous series and open to some interpretation.

2.1 GPlates

GPlates (https://www.gplates.org) is free software that reconstructs the movement of tectonic plates through geological history. Using models of plate movement, the configuration can be reconstructed at arbitrary points in the past. Onto this, the palaeoenvironments of the formations at their current positions can be overlain, and so the positions of landmasses, mountains, coastlines, and marine settings can be plotted long back in time.

The download from GPlates provides lots of example data and projects to show the various plates, tectonics, plumes and other information that can be included. New data are being added as research is published. Try, for instance, the DataBundleForNovices project (GplatesSampleData/DataBundleForNovices) to see a combined model of modern continents and ocean spreading over the last 410 million years.

The main GPlates window showing the DataBundleForNovices loaded.

Layers can be selected in this Gplates window to show and hide different features.

There is also a GPlates Web Service (GWS; http://gws.gplates.org) that can provide some of the same data without the software. I’ll look briefly at this first.

2.2 Exercises 1

1. Find an example of a palaeogeographical map for your favourite time period.
2. Have a look in the folder GplatesSampleData to see what other projects are included.

3 Automatic map plotting in R

A couple of packages already exist for using GPlates reconstructions to create palaeogeographical maps in R:

• NonaR/paleoMap (https://github.com/NonaR/paleoMap)
  • NB Has not been updated in several years.
• LunaSare/gplatesr (https://github.com/LunaSare/gplatesr)
  • More recently updated, and similar to above.

Both of these download from GWS. I’ve borrowed some of their functions for this lab group.

If you’re unfamiliar with RStudio or RMarkdown (.Rmd files), then the code sections are placed in ‘chunks’ like the one below. You can run individual lines by Cmd+Enter or Ctrl+Enter, or run the whole chunk using the little green arrow to the top right of each chunk. Try with this chunk below if you want.

library(broom)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(forcats)
library(ggplot2)

Use suppressPackageStartupMessages() to eliminate package startup
messages

library(ggthemes)
library(mapproj)

Loading required package: maps

library(purrr)

Attaching package: 'purrr'

The following object is masked from 'package:maps':

    map

library(readr)
library(rgdal)

Loading required package: sp

rgdal: version: 1.5-23, (SVN revision 1121)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
Linking to sp version:1.4-5
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

library(stringr)
library(tibble)
# library(tidyverse) # a quick way to load packages, but slow for Binder.

list.files("functions/", pattern = "\\.R", full.names = TRUE) %>%
  walk(source)

NB not all of the chunks below should be run on Binder as some of the code will take a long time or fail to execute.

3.1 GPlates Web Service

GWS can reconstruct the locations of points, coastlines, plates, paths, and features back into the past (https://github.com/GPlates/gplates_web_service_doc/wiki). It can also use a variety of different reconstructions. Most relevant is perhaps the Golonka (2007) model, which is often used for palaeocoordinate reconstruction in the Palaeobiology Database (PBDB).

The process is:

1. Download outlines of the continental coastlines.
2. Download outlines of the tectonic plates.
3. Read and convert these into a format that R and ggplot2 can use,

• this uses the packages rgdal, readr, and broom.

4. Plot,

• in this case with gpplot.

Example code only. Do not run. The code below does these steps by downloading the data, which will likely take some time.

#### DO NOT RUN THIS CODE ####
coastline_gws_url <-
  "http://gws.gplates.org/reconstruct/coastlines/?time=155&model=GOLONKA"
polygons_gws_url <-
  "http://gws.gplates.org/reconstruct/static_polygons/?time=155&model=GOLONKA"

kimmeridgian_coastlines <-
  rgdal::readOGR(coastline_gws_url) %>%
    broom::tidy()
kimmeridgian_polygons <-
  rgdal::readOGR(polygons_gws_url) %>%
    broom::tidy()

This chunk below replicates the above using data I’ve already downloaded from GWS and stored in the GitHub repository, and can be used on Binder. It has the same steps of loading and tidying the data, then includes the commands to plot with ggplot2.

Feel free to run this code.

# Load and tidy the data 
kimmeridgian_coastlines <-
  readOGR(
    "data/GWS/Matthews_etal_GPG_2016_Coastlines_reconstructed_155.00Ma.gmt"
  ) %>%
  tidy()
kimmeridgian_polygons <-
  readOGR(
    "data/GWS/Matthews_etal_GPG_2016_Polygons_reconstructed_155.00Ma.gmt"
  ) %>%
  tidy()

ggplot() +
  geom_map(
    data = kimmeridgian_polygons,
    map = kimmeridgian_polygons,
    aes(x = long, y = lat, map_id = id),
    fill = "#D8D8D8"
  ) +
  geom_map(
    data = kimmeridgian_coastlines,
    map = kimmeridgian_coastlines,
    aes(x = long, y = lat, map_id = id),
    colour = "#222222", fill = NA, size = 0.3
  ) +
  coord_map("mollweide") +
  map_border() +
  theme_map()

Outlines of continental plates in the Kimmeridgian (155 Ma). Plates in the continental model are shaded grey, while the modern coastlines are outlined. Data downloaded from the GPlates web service.

When plotting, the data are added as layers, each plotting over the other: plate polygons at the bottom then coastlines on top. I’ve also set the map projection to be Mollweide, but you can play with this and change it to others included in the mapproj package function mapproject. Options include "Gilbert", "fisheye", and "hex". (Use, for example, coord_map("fisheye", n = 0.7) to add extra arguments.)

Additionally I’ve added a border (map_border(), in the functions folder) and used a basic theme (theme_map()) to get a plain background.

3.2 How R plots these map data

The positions and extent of map areas from data in GPlates are stored as a series of polygons. Each region is outlined by a series of points that form a closed shape, filled with a particular colour. This is exported from GPlates as OGR-GMT format, with file extension .gmt, which lists each point that marks the corner of a polygon. In R, broom converts this to a tibble (similar to a data.frame) that has the x and y locations of each point, and a series of columns that identify which polygon that point belongs to. These ID columns are important to make sure that polygons contain the correct points and are closed properly. Without these, there would only be a single shape that would stretch across the whole globe and not make much sense.

3.3 Not quite right

GWS is convenient for data and access, but the results in the plot above show the current continental coastlines and the plates. This does not reflect the ancient coastlines that were really after. Also, many of the countries are cross-cut by the polygons that form the sections in the model, which are untidy and typically not visible.

We need to go to the next level: true palaeogeographical reconstructions. Fortunately GPlates offers that too.

3.4 Exercises 2

1. Have a look at the GWS documentation. What address would you use to get reconstructed coastlines in the early Carboniferous using the Matthews et al. (2016) model?
2. Open an example of the OGR-GMT data from the GWS, either from the data/GWS folder or by downloading it. Can you see how the polygons are stored? How are different polygons identified?
3. Now look at the tidied version in R, i.e. kimmeridgian_coastlines or kimmeridgian_polygons. By default this will only show the first 10 rows, but you can see more with, for example, head(kimmeridgian_polygons, n = 50L). Can you see how the OGR-GMT data is transferred to this tidy format?
4. Play around with the inputs to the plot above. How messy or colourful can you make this map? What options are there in theme_map()? (Hint: use ?theme_map to see.)

4 ‘True’ Palaeogeographical Outlines

Among the projects included with the GPlates download is the set of palaeogeographical reconstructions by Cao et al. (2017) for the last 402 my. These are based, and so should align nicely, with the Golonka (2007) model that is used for data in the PBDB, which we’ll get to below. These data are in GPlatesSampleData/FeatureCollections/Palaeogeography/Global/.

Cao et al. (2017) used a combination of previous palaeogeographical reconstructions, and updated the coastlines, mountains, and extent of land and ice caps using occurrence data from the PBDB. In a few locations, what had been reconstructed as land produced marine fossils, when accounting for the movement of plates, and vice versa.

These reconstructions include four layers that have the separated palaeoevonvironmental regions:

• Shallow marine: shelf seas.
• Land: terrestrial habitats outlined by coastlines.
• Mountain: areas of mountain building and uplift.
• Ice cap: polar regions that have substantial evidence for long term presence of ice.

These are given the abbreviations ‘sm,’ ‘l,’ ‘m’ and ‘i.’ It’s worth noting that the way these are built means that the layers have to be plotted in the order above – shallow marine first, ice cap last. There are various areas that overlap (Australia especially) so the order of overlay becomes important in making the correct boundaries visible.

4.1 Export from GPlates

The Cao et al. (2017) reconstructions aren’t available from GWS, but can be exported from GPlates directly using the menu item Reconstruction > Export… Like above I’ve already done this for an example data set, in this case using the Late Jurassic (Kimmeridgian, 155 Ma) as an example. To make your own, once in the export window:

1. Select the time of export (can be set from the main window time).
   • Also an option to export a series at regular intervals.
2. The layers to export
   1. select Add Export… under Export Data
   2. choose Reconstructed Geometries
   3. output as OGR-GMT
   4. choose single or multiple files
      • I prefer multiple: you get one file for each of sea/land/mountain/ice. This is the form that works with the code below.
   5. OK
3. Change Target directory: to somewhere you can find.
4. Begin Export

Export options from GPlates.

Layer selection in GPlates export. Select Reconstructed geometries, OGR-GMT, and multiple layers.

The example data are in data/palaeogeography as three .gmt files that describe the polygon outlines for each of the layers. Note that as there are no ice caps at this time, there are only three files/layers. Look at the Carboniferous of later Palaeogene to get some ice caps.

4.2 Read in the map data

I find it useful to tie together the layers and colours from the beginning of importing – this ensures keeping the layers in the order for plotting. The data files are:

• lm_402_2_reconstructed_155.00Ma.gmt: landmass
• m_402_2_reconstructed_155.00Ma.gmt: mountain
• sm_402_2_reconstructed_155.00Ma.gmt: shallow marine

This short code snippet sets lists the layer names and assigns the colours for each layer on the map – they are similar to those shown in GPlates, but not quite the same. I’ve also used #DAD3FF for ice caps when needed. (Technically this isn’t used for the colours, that’s in the palaeogeog_map_niceties function, this is somewhat a remnant of a previous way that I made this.)

map_layers <-
  c(
    "Land"           = "#FFD23A",
    "Mountain"       = "#FF8D51",
    "Shallow marine" = "#45D8FF"
  )

With that ready, the following code reads and tidies the data (as for above), then adds a column (geog_layer) to identify the data layer (shallow marine, land, mountain) and joins all the data together. It’s very tidyverse based, so I’d suggest looking at https://www.tidyverse.org to see the rationale and how they link using the pipe (%>%). You can also see more about individual functions in the R help.

The final step is to create a column (polygon_id) to identify each region on the, which come after mountain map. I had a long issue where different continents would join up, or Australia’s land would always appear plotted underneath the sea. As described above, this is because the shapes on the map are plotted as separate polygons, and the IDs are repeated between each of the layers. We’ll use the polygon_id column to make sure that each polygon has a unique ID by joining the layer, group and id name from the tidied data.

polygon_data <-
  list.files("data/palaeogeography/", pattern = ".gmt", full.names = TRUE) %>%
    map(readOGR) %>%
    map(tidy) %>%
    map2(
      names(map_layers),
      ~ add_column(.x, geog_layer = .y)
    ) %>%
    bind_rows() %>%
    mutate(
      polygon_id = str_c(id, group, geog_layer) %>% as_factor(),
      geog_layer = factor(geog_layer, labels = names(map_layers))
    )

OGR data source with driver: OGR_GMT 
Source: "/Users/glbcm/GitHub/palaeomap_example/data/palaeogeography/lm_402_2_reconstructed_155.00Ma.gmt", layer: "lm_402_2_reconstructed_155.00Ma"
with 361 features
It has 20 fields
OGR data source with driver: OGR_GMT 
Source: "/Users/glbcm/GitHub/palaeomap_example/data/palaeogeography/m_402_2_reconstructed_155.00Ma.gmt", layer: "m_402_2_reconstructed_155.00Ma"
with 195 features
It has 20 fields
OGR data source with driver: OGR_GMT 
Source: "/Users/glbcm/GitHub/palaeomap_example/data/palaeogeography/sm_402_2_reconstructed_155.00Ma.gmt", layer: "sm_402_2_reconstructed_155.00Ma"
with 808 features
It has 20 fields

Regions defined for each Polygons
Regions defined for each Polygons
Regions defined for each Polygons

4.3 Arranging the map layers

At this point it’s important to have the map layers in the correct order, so that marine is plotted under land, which is plotted under mountain.

This piece of code does that for you, ensuring shallow marine is first (it was third in map_layers above) then assigning this order to the polygon_id column. It converts the polygon_id column to a factor and re-orders thee levels based on the reordered map_layers.

id_level_order <-
  map_layers[c(3, 1, 2)] %>%
    map(
      function(layr) polygon_data$polygon_id %>% levels() %>% str_which(layr)
    ) %>%
    unlist()

polygon_data <-
  polygon_data %>%
    mutate(
      polygon_id = fct_relevel(polygon_id, levels(polygon_id)[id_level_order])
    )

4.4 Plotting the map

Unfortunately, using the same plot layer method as for the GWS data above is a little excessive. You would have to add a new annotation for each layer (just about okay for three, tedious for more), which then multiplies when plotting reconstructions of many times.

# DO NOT RUN THIS CODE
ggplot() +
  geom_map() +
  geom_map() +
  geom_map() + …

Instead, you can use the function geom_polygon and have the different layers plotted automatically with the colours assigned above. The fill is the most important argument as that assigns the main bulk of the colour. The colour argument is the outline, which you shouldn’t use bas that will outline the map layers and change their shape. The colour argument will be used when plotting occurrence points (see below). Also having the ordered and correctly grouped polygon_id is important here as that determines how the individual polygons are grouped and closed, preventing shapes from spanning the globe, or incorrectly overlapping.

I’ll point out here that there’s quite a lot going on in the palaeogeog_map_niceties function at the end. This is a collection of theme options and a colour scale to show the map properly. In particular this adds the outline and the lines of the equator and tropics. I guess one thing to note is that the ‘tropics,’ in a biodiversity sense, may not be the same in the past as the present. Increased temperatures can change the latitudinal biodiversity gradient. Also, obliquity and precession change the latitudes of the astronomical tropics through time.

# dev.new(width = 7, height = 4)
map_plot <-
  ggplot() +
    geom_polygon(
      data = polygon_data,
      aes(
        x      = long,
        y      = lat,
        fill   = geog_layer,
        group  = polygon_id
      ),
      colour = NA
    ) +
    coord_map("mollweide") +
    theme_map() +
    palaeogeog_map_niceties()
map_plot

Palaeogeography in the Kimmeridgian (155 Ma). Palaeogeographical reconstruction from Cao et al. (2017). Lines of latitude indicate the locations of the Equator and Tropics (±23.5°).

You can try changing the look by modifying the projection (coord_map) and play about with some of the colours and markings by editing the function palaeogeog_map_niceties (functions/palaeogeog_map_niceties.R). You can also add further annotations (https://ggplot2-book.org/annotations.html#annotations) or change labels as you wish.

4.5 Going further 1: locating with modern coastlines

You may also want to add outlines of modern coastlines to show where these are and help locating. Again these are not available from GWS, but can be exported from GPlates. In GPlates, from File > Open Feature Collection… go to GplatesSampleData/FeatureCollections/Coastlines and load Matthews_etal_GPC_2016_Coastlines_Polyline.gpmlz. This adds a new layer into GPlates with areas of modern coastline in their relative positions for the time. Not all the modern coastlines are there as not all of these have certain locations, or rock to be found at – these come and go at different time points.

Like above, you can export this layer from GPlates. In this repo, I’ve already done that to data/coastlines/, which you can load with this code.

modern_coastlines <-
  readOGR("data/coastlines/Matthews_etal_GPC_2016_Coastlines_Polyline_reconstructed_155.00Ma.gmt") %>%
    tidy() %>%
    add_column(geog_layer = "Modern coastlines") %>%
    mutate(polygon_id = str_c(id, group, geog_layer) %>% as_factor())

This layer needs the same organising of IDs as the palaeogeography above, but once done can be easily added to map_plot above. This is where it’s useful to save your ggplot to an object alongside printing it – this you can just add new layers on top of the base map.

map_plot +
  geom_path(
    data = modern_coastlines,
    aes(
      x = long,
      y = lat,
      group = polygon_id,
    ),
    colour = "#888888", size = 0.3
  )

Palaeogeography in the Kimmeridgian (155 Ma). Outlines of modern coastlines (where known) are included in grey. Palaeogeographical reconstruction from Cao et al. (2017). Lines of latitude indicate the locations of the Equator and Tropics (±23.5°).

4.6 Exercises 3

1. Try plotting this map with different projections (see above and ?mapproject). Does this work well for all projections? Why?
2. Look at the code behind palaeogeog_map_niceties() (in functions/palaeogeog_map_niceties.R). Modify this to add different lines of latitude. Can you add lines of longitude too? What about changing the colour?

• ggplot objects can have new layers added to them one-by-one by using +, but you can also combine them into function containing a list of layers that can be added in a single line.
• Hint: to change the colour of the map layers, use fill.

3. Can you try creating your own palaeogeographical map for your favourite time? This should be done on your own computer. Follow the instructions above to export data from GPlates then plot in R, making sure you change all the file and object names.

5 Adding Fossil Occurrences

5.1 Getting occurrence data

Now that you’ve spent lots of time making a pretty map, let’s plot some fossils on top of it.

Of course, being me and having used Jurassic data, we’re going to plot some ichthyosaurs. In this instance, the occurrence data comes from the PBDB, downloading ichthyosaur occurrences with these filters:

• Callovian–Tithonian (166-145 Ma)
• all levels of taxonomy – no filtering.

The PBDB has its own web API (https://paleobiodb.org/data1.2) from which you can get occurrences, collections, taxonomy and other things. In the code chunk below I’ve crafted the URL to download the data, but this is already in the repo so doesn’t need to be run.

Do not run the following chunk, the data is already in the repo.

#### DO NOT RUN THIS CHUNK ####
pbdb_url <-
  "https://paleobiodb.org/data1.2/occs/list.csv?base_name=Ichthyosauromorpha&interval=Callovian,Oxfordian,Kimmeridgian,Tithonian&show=paleoloc"

occ_ichthyosaurs <-
  read_csv(pbdb_url)

The part of the URL show=paleoloc is the important bit as the PBDB has already computed the palaeocoordinates of most occurrences using their modern location and the ages of the occurrences. Conveniently using the same default GPlates model as the palaeogeography.

5.2 Plotting occurrences

With that all done, use the chunk below to add the occurrences onto the base map. These are plotted with geom_point(), which you may be familiar with, setting the x and y values to the palaeolongitude and palaeolatitude respectively.

occ_ichthyosaurs <-
  read_csv("data/occurrences/2021-02-18-ichthyosaur-occurrences.csv")

── Column specification ────────────────────────────────────────────────────
cols(
  .default = col_character(),
  occurrence_no = col_double(),
  reid_no = col_double(),
  flags = col_logical(),
  collection_no = col_double(),
  identified_no = col_double(),
  accepted_no = col_double(),
  max_ma = col_double(),
  min_ma = col_double(),
  reference_no = col_double(),
  paleolng = col_double(),
  paleolat = col_double()
)
ℹ Use `spec()` for the full column specifications.

occ_plot <-
  map_plot +
    geom_point(
      data = occ_ichthyosaurs,
      aes(
        x = paleolng,
        y = paleolat
      )
    )
occ_plot

Occurrences of ichthyosaurs from the Callovian–Tithonian. Palaeogeographical map shows the distribution of land in the Kimmeridgian (155 Ma) reconstructed by Cao et al. (2017). Occurrences of Ichthyosauromorpha are from the Palaeobiology Database.

The more eagle-eyed of you may notice that some of these marine ichthyosaur occurrences are located resolutely in the middle of landmasses. There are a few possible reasons for that:

1. The palaeogeographical reconstruction is wrong.
2. The modern location for the occurrence is wrong, so the palaeocoordinates are wrong.
3. The age of the occurrence is wrong, so its reconstructed position is wrong.

These three are worth considering further.

Incorrect palaeogeographical reconstruction is perhaps the most obvious answer: these reconstructions are difficult and there will always be some uncertainty with reconstructing the earth 155 million years ago. However, the method of Cao et al. (2017) included palaeoenvironments reconstructed using PBDB data to improve confidence in the borders of sea and land. The palaeogeography should match the PBDB well then, particularly in well-studied and sampled periods and locations. There are, however, some parts that cannot certainly reconstructed because of a lack of data, rock availability, or accessibility.

Incorrect modern location may be the easiest to discount as localities tend to be spatially constrained and identifiable, even in older literature. This being the case, there should only be minor error in the locations of occurrences (~10s km).

Incorrect age assignment is probably the most likely cause, or perhaps better imprecise ages. The palaeocoordinates calculated for each occurrence in the PBDB use the midpoint of its temporal range. Sometimes this range can be quite large, such as the several species of Undorosaurus that cover the whole Late Jurassic (163.5–145 Ma in GTS 2020/03). Those specimens from Patagonia that are plotted in the middle of the land are also similarly imprecisely dated. Their palaeocoordinates are reconstructed to the middle Kimmeridgian, whereas they are likely later (given the occurrence of ichthyosaurs from similar locations) when there was a sea there to swim in.

As the ichthyosaur occurrences included over over 20 my of earth history there will be changes in palaeogeography, and some potentially imprecise dating.

The ways around this are to have more time slices with the occurrences plotted, but this complicates the code and any plots, and such precise reconstructions may not be available. Alternatively, aim to use only the most precisely dated occurrences, if possible, but this can remove lots of useful data. Ultimately, you may just have to accept that some marine occurrences may appear on land, or vice versa.

5.3 Going further 2: separating occurrences with facets

The final thing that you can try is the faceting power of ggplot2 to separate the different groups of occurrences. Below, this code splits out the occurrences by their accepted taxonomic rank. Most occurrences are just non-diagnostic vertebral material, so can’t be assigned to a family, genus or species. These are listed under ‘unranked_clade.’

Conveniently, creating facets is as simple as adding a single line.

occ_plot +
  facet_wrap(vars(accepted_rank), ncol = 2)

Occurrences of ichthyosaurs in the Callovian–Tithonian separated by identified rank. Palaeogeographical maps shows distribution of land in the Kimmeridgian (155 Ma) reconstructed by Cao et al. (2017).

5.4 Exercises 4

1. Look at the PBDB web API (https://paleobiodb.org/data1.2). What other options and data sets can you access with this? Can you craft a URL that will download occurrence data of Tithonian bivalves.

• This download may not work on Binder, so you should do it in R on your own computer.

2. Try adding colour or shapes to you occurrence points using the options in geom_point.
3. Download and plot the palaeogeographical map and occurrences for your favourite group in your favourite stage. Again this should be done on your own computer. Beware it may take a while if you want to plot many occurrences.
4. It’s also useful to separate out occurrences on the map, for instance, can you colour occurrences within the outside the tropics differently? say red in the tropics (±23.5°) and blue outside.

• Hint: you can use the palaeocoordinates of the PBDB data to assign a group to each occurrence, use this to colour the points in geom_point.
• Hint: use scale_colour_discrete to colour the points. See the code of palaeogeog_map_niceties (functions/palaeogeog_map_niceties.R) for an example of doing this.

6 Suggested solutions to exercises

Not all of the exercise questions have a single obvious solution, particularly ones where you pick your favourite period. But here are some code suggestions for you to check against.

6.1 Exercises 1

1. Find an example of a palaeogeographical map for your favourite time period.

• You can have a look at Deep Time Maps or Paleomap for easy examples, or you can find others in papers.

2. Have a look in the folder GplatesSampleData to see what other projects are included.

• The GplatesSampleData/FeatureCollections directory is the one to look at. It includes:
  • Other plate reconstructions.
  • Ocean boundaries.
  • Plate flow lines.
  • Large igneous provinces … among other things.

6.2 Exercises 2

1. Have a look at the GWS documentation. What address would you use to get reconstructed coastlines in the early Carboniferous using the Matthews et al. (2016) model?

• "https://gws.gplates.org/reconstruct/coastlines/?&time=350&model=MATTHEWS2016_mantle_ref"
• You could also use "MATTHEWS2016_pmag_ref". These are two different reconstruction methods.

2. Open an example of the OGR-GMT data from the GWS, either from the data/GWS folder or by downloading it. Can you see how the polygons are stored? How are different polygons identified?

• The metadata at the top has lines beginning with # @, including describing how the data for each polygon is organised.
• Each polygon has a line of metadata describing it and @P indicating it’s a polygon.
• Sometimes without metadata for each polygon, they are indicated with numbers: @P1, @P2 etc.

3. Now look at the tidied version in R, i.e. kimmeridgian_coastlines or kimmeridgian_polygons. By default this will only show the first 10 rows, but you can see more with, for example, head(kimmeridgian_polygons, n = 50L). Can you see how the OGR-GMT data is transferred to this tidy format?

• The data is stored in a tibble with each point in a row and columns for the longitude and latitude of the points.
• The order that the points should be plot and the piece, group, and id of the polygon that they belong to are also included to separate the shapes.
• These ID columns are taken from the metadata for each polygon define in the input OGR-GMT file.

4. Play around with the inputs to the plot above. How messy or colourful can you make this map? What options are there in theme_map()? (Hint: use ?theme_map to see.)

• I’ll leave this to your imagination.

6.3 Exercises 3

1. Try plotting this map with different projections (see above and ?mapproject). Does this work well for all projections? Why?

• Some of the projections move the line of 0° longitude away from the centre or have polygons that cross from one side of the map to the other. This can sometimes lead to shapes being drawn right across the map. I’m not sure whether the best way to counteract this is in redrawing the polygons or something to do with the map projection.

2. Look at the code behind palaeogeog_map_niceties() (in functions/palaeogeog_map_niceties.R). Modify this to add different lines of latitude. Can you add lines of longitude too? What about changing the colour?

• In palaeogeog_map_niceties I stored the lines of latitude that I want as a tibble that describe line segments. New lines of latitude can be added by creating a new row with the value y set to your latitude. The lines are described as going from -180°–+180° longitude at a specific single latitude.
• For lines of longitude, I’d suggest using a new geom_segment layer. Modify the code so that the y value covers the range of latitudes –90°–+90° at a specific longitude. Here’s an example drawing lines of longitude at -50°, 0° and 72°:

geom_segment(
  data = tibble(
    x = c(-50, 0, 72),
    y = rep(-90, 3),
    yend = rep(90, 3)
  ),
  aes(x = x, xend = x, y = y, yend = yend),
  colour = "grey70", size = 0.3
)

3. Can you try creating your own palaeogeographical map for your favourite time? This should be done on your own computer. Follow the instructions above to export data from GPlates then plot in R, making sure you change all the file and object names.

• here you’ll have to modify the code above for your own GPlates data. Remember to add in the ice caps to the layer if you pick a time with these, such as the Carboniferous or the Neogene.

6.4 Exercises 4

1. Look at the PBDB web API (https://paleobiodb.org/data1.2). What other options and data sets can you access with this? Can you craft a URL that will download occurrence data of Tithonian bivalves.

• “https://paleobiodb.org/data1.2/occs/list.csv?base_name=Bivalvia&interval=Tithonian&show=paleoloc”

2. Try adding colour or shapes to your occurrence points using the options in geom_point.

• Use the colour and shape arguments.

3. Download and plot the palaeogeographical map and occurrences for your favourite group in your favourite stage. Again this should be done on your own computer. Beware it may take a while if you want to plot many occurrences.

• Use the code above, but change the relevant URLs and data files to your own GPlates and PBDB data. Don’t forget to add ice caps if you have them.
• Hopefully the PBDB API will be working as it wasn’t when I was writing this!

4. It’s also useful to separate out occurrences on the map, for instance, can you colour occurrences within the outside the tropics differently? say red in the tropics (±23.5°) and blue outside.

• You’ll need to assign a column of identifiers for tropical and non-tropical occurrences. Below is an example using the dplyr function case_when.
• With that done it’ll fit right into the plotting code from earlier, but you will need to specify the colour argument. If you want to make it blue and red then you can add a scale.

ex_occ_ichthyosaurs <-
  occ_ichthyosaurs %>%
    mutate(
      tropical = case_when(
        paleolat <= 23.5 & paleolat >= -23.5 ~ "tropical",
        paleolat > 23.5 | paleolat < -23.5   ~ "non-tropical"
      ) %>%
        as_factor()
    )

ex_occ_plot <-
  map_plot +
    geom_point(
      data = ex_occ_ichthyosaurs,
      aes(
        x = paleolng,
        y = paleolat,
        colour = tropical
      ),
      na.rm = TRUE
    ) +
    scale_discrete_manual(
      values = 
        c(
          "tropical" = "red",
          "non-tropical" = "blue"
        ),
        aesthetics = c("colour"),
        name = "Latitude"
      )
ex_occ_plot

Occurrences of ichthyosaurs in the Callovian–Tithonian indicated by latitude region. Palaeogeographical maps shows distribution of land in the Kimmeridgian (155 Ma) reconstructed by Cao et al. (2017).

7 References

Bradshaw, M J, John Christopher Wolverson Cope, D W Cripps, D T Donovan, Michael Kingsley Howarth, P F Rawson, I M West, and William Albert Wimbledon. 1992. “Jurassic.” Geological Society, London, Memoirs 13 (1): 107–29. https://doi.org/10.1144/GSL.MEM.1992.013.01.12.

Cao, Wenchao, Sabin Zahirovic, Nicolas Flament, Simon Williams, Jan Golonka, and R. Dietmar Müller. 2017. “Improving global paleogeography since the late Paleozoic using paleobiology.” Biogeosciences 14 (23): 5425–39. https://doi.org/10.5194/bg-14-5425-2017.

Gao, Ting, Da-Qing Li, Long-Feng Li, and Jing-Tao Yang. 2019. “The First Record of Freshwater Plesiosaurian from the Middle Jurassic of Gansu, NW China, with Its Implications to the Local Palaeobiogeography.” Journal of Palaeogeography, August. https://doi.org/10.1186/s42501-019-0043-5.

Golonka, J. 2007. “Phanerozoic Paleoenvironment and Paleolithofacies Maps: Mesozoic.” Geologia 33 (2): 211–64. https://www.infona.pl//resource/bwmeta1.element.baztech-article-AGH5-0014-0139.

Matthews, Kara J., Kayla T. Maloney, Sabin Zahirovic, Simon E. Williams, Maria Seton, and R. Dietmar Müller. 2016. “Global Plate Boundary Evolution and Kinematics Since the Late Paleozoic.” Global and Planetary Change 146 (November): 226–50. https://doi.org/10.1016/j.gloplacha.2016.10.002.

Zverkov, Nikolay G., Dmitry V. Grigoriev, and Igor G. Danilov. 2020. “Early Jurassic Palaeopolar Marine Reptiles of Siberia.” Geological Magazine, 1–18. https://doi.org/10.1017/S0016756820001351.